home *** CD-ROM | disk | FTP | other *** search
/ Gold Medal Software 2 / Gold Medal Software Volume 2 (Gold Medal) (1994).iso / windows / win31 / macsyma.arj / MACSDEMO.EXE / LINDST.OUT < prev    next >
Text File  |  1993-09-14  |  8KB  |  119 lines
  1.  
  2. (c1) /*
  3.                           DEMONSTRATION OF PACKAGE "LINDST"
  4.  
  5.                                                                               
  6.    This file shows how to use the LINDSTEDT program to compute periodic 
  7.    solutions of weak polynomial perturbations a single harmonic oscillator.
  8.    For a more detailed demonstration, see the demo file LINDST1.  
  9.    If the perturbation depends explicitly on the independent variable then
  10.    it must be a periodic dependence.
  11.  
  12.    The program attempts to compute all periodic solutions of a planar system
  13.    of the form
  14.  
  15.                P*X' =  A*X + B*Y + F0(T) + E*F1(T,X,Y) + E^2*F2(T,X,Y) + ...
  16.                Q*Y' =  C*X - D*Y + G0(T) + E*G1(T,X,Y) + E^2*G2(T,X,Y) + ...
  17.  
  18.    where the coefficients A,B,C,D,P,Q are such that the origin is an elliptic 
  19.    fixed  point of the equation when E=0 and F0=G0=0.  It will also accept a 
  20.    single second-order equation of the form
  21.  
  22.               A*X'' + B*X + F0(T) + E*F1(T,X,X') + E^2*F2(T,X,X') + ... = 0
  23.  
  24.    provided that B/A is positive if it has a numerical value.
  25.  
  26.    A first-order system must be specified as a list, and a single second-order
  27.    equation must NOT be given as a list.  The perturbation parameter and 
  28.    truncation order must also be given.  An optional argument specifying the
  29.    initial conditions which will be satisfied at zeroth order can be given.
  30. */
  31. if get('lindstedt,'version)=false then load(lindst)$
  32. C:\MACSD2B\share\LINDST.fas being loaded.
  33.  
  34. (c2) /*
  35.    If the equations are independent of the independent variable, denoted by T
  36.    in the examples given above, then the program attempts to compute all 
  37.    unique periodic solutions which have non-negative zeroth order amplitude.
  38.    (Setting the option variable LINDSTEDT_DELETE_LC [TRUE] to FALSE forces 
  39.    the program to return all solutions.)
  40.  
  41.    The initial condition (A,B) at T=0 can be given as an optional argument.  
  42.    This can be used to compute a particular periodic solution when the 
  43.    equations have a continuous family of periodic solutions.  It should not 
  44.    be used when attempting to compute a limit cycle, since the specification 
  45.    of initial conditions will be inconsistent in this case.
  46.  
  47.    If the perturbation depends explicitly on the independent variable, the 
  48.    program will attempt to compute all periodic solutions.  However, it is 
  49.    often the case that the resulting secular equations are too complex for 
  50.    MACSYMA to solve exactly.
  51. */
  52.  
  53. disp(" ")$
  54. |$label(0,15,Times New Roman,) 
  55.  
  56.  
  57. (c3) /* 
  58.    The first example is the familiar Van der Pol oscillator, which is known to
  59.    have a limit cycle for all nonzero E.  The limit cycle will be computed for
  60.    the original second-order equation, and again for the equivalent first-
  61.    order system. */
  62.  
  63. (kill(x,y,t,e,vdp,vdp_system,duffing),depends([x,y],t),
  64. vdp: 'diff(x,t,2)+x-e*'diff(x,t)*(1-x^2)=0);
  65. |$label(0,15,Times New Roman,$(d3$))$q($sup(d,2)x,d$sup(t,2))$hinge()$in( - )e$in( )$paren(1$in( - )$sup(x,2),$(,$))$in( )$q(dx,dt)$hinge()$in( + )x$hinge()$in( = )0
  66.  
  67. (c4) /* Now compute the limit cycle to O(E^2). This will take a minute...*/
  68.  
  69. lindstedt(vdp,e,2);
  70. C:\MACSD2B\library1\algsys.fas being loaded.
  71. C:\MACSD2B\library1\grobner.fas being loaded.
  72. C:\MACSD2B\library1\result.fas being loaded.
  73. |$label(0,15,Times New Roman,$(d4$))$open([)$open([)$open([)$in( - )$q($paren(5$in( )cos$paren(5$in( )%tau)$in( - )18$in( )cos$paren(3$in( )%tau)$in( + )12$in( )cos$paren(%tau),$(,$))$in( )$sup(e,2),96)$hinge()$in( - )$q($paren(sin$paren(3$in( )%tau)$in( - )3$in( )sin$paren(%tau),$(,$))$in( )e,4)$hinge()$in( + )2$in( )cos$paren(%tau)$close(])$ina($, )$hinge()%tau$hinge()$in( = )$open($()1$hinge()$in( - )$q($sup(e,2),16)$close($))$hinge()$in( )t$close(])$close(])
  74.  
  75. (c5) /* 
  76.    The program returns a list of solutions.  Each solution consists of the 
  77.    periodic orbit, the time scaling, and any unresolved secular equations,
  78.    if they are present. In this case, there are no unresolved secular 
  79.    equations. 
  80.  
  81.    We next apply LINDSTEDT to a first-order system which is equivalent to the
  82.    Van der Pol equation. */
  83.  
  84. vdp_system: ['diff(x,t)=y,'diff(y,t)=-x+e*y*(1-x^2)];
  85. |$label(0,15,Times New Roman,$(d5$))$open([)$q(dx,dt)$hinge()$in( = )y$ina($, )$hinge()$q(dy,dt)$hinge()$in( = )e$in( )$paren(1$in( - )$sup(x,2),$(,$))$in( )y$hinge()$in( - )x$close(])
  86.  
  87. (c6) lindstedt(vdp_system,e,2);
  88. C:\MACSD2B\share\blinalgl.fas being loaded.
  89. |$label(0,15,Times New Roman,$(d6$))$open([)$open([)$open([)$in( - )$q($paren(5$in( )cos$paren(5$in( )%tau)$in( - )18$in( )cos$paren(3$in( )%tau)$in( + )12$in( )cos$paren(%tau),$(,$))$in( )$sup(e,2),96)$hinge()$in( - )$q($paren(sin$paren(3$in( )%tau)$in( - )3$in( )sin$paren(%tau),$(,$))$in( )e,4)$hinge()$in( + )2$in( )cos$paren(%tau)$ina($, )$q($paren(25$in( )sin$paren(5$in( )%tau)$in( - )54$in( )sin$paren(3$in( )%tau)$in( + )24$in( )sin$paren(%tau),$(,$))$in( )$sup(e,2),96)$hinge()$in( - )$q($paren(3$in( )cos$paren(3$in( )%tau)$in( - )3$in( )cos$paren(%tau),$(,$))$in( )e,4)$hinge()$in( - )2$in( )sin$paren(%tau)$close(])$ina($, )$hinge()%tau$hinge()$in( = )$open($()1$hinge()$in( - )$q($sup(e,2),16)$close($))$hinge()$in( )t$close(])$close(])
  90.  
  91. (c7) /* 
  92.    Since the equation was given as a system, the program returns both the X
  93.    and Y components of the solution. */
  94.  
  95. disp(" ")$
  96. |$label(0,15,Times New Roman,) 
  97.  
  98.  
  99. (c8) /* 
  100.    Next, consider the Duffing equation.  It is known that all solutions of 
  101.    this equation are periodic and that they fill the phase plane.  We will 
  102.    compute the periodic solution through O(E^2) with initial conditions (A,B).
  103. */
  104.  
  105. duffing: 'diff(x,t,2)+x+e*x^3=0;
  106. |$label(0,15,Times New Roman,$(d8$))$q($sup(d,2)x,d$sup(t,2))$hinge()$in( + )e$in( )$sup(x,3)$hinge()$in( + )x$hinge()$in( = )0
  107.  
  108. (c9) lindstedt(duffing,e,2,[a,b]);
  109. |$label(0,15,Times New Roman,$(d9$))$open([)$open([)$open([)$q($paren($paren(sin$paren(5$in( )%tau)$in( + )48$in( )sin$paren(3$in( )%tau)$in( + )271$in( )sin$paren(%tau),$(,$))$in( )$sup(b,5)$in( + )$paren(5$in( )cos$paren(5$in( )%tau)$in( + )108$in( )cos$paren(3$in( )%tau)$in( - )113$in( )cos$paren(%tau),$(,$))$in( )a$in( )$sup(b,4)$in( + )$paren($in( - )10$in( )sin$paren(5$in( )%tau)$in( + )12$in( )sin$paren(3$in( )%tau)$in( + )854$in( )sin$paren(%tau),$(,$))$in( )$sup(a,2)$in( )$sup(b,3)$in( + )$paren($in( - )10$in( )cos$paren(5$in( )%tau)$in( + )180$in( )cos$paren(3$in( )%tau)$in( - )170$in( )cos$paren(%tau),$(,$))$in( )$sup(a,3)$in( )$sup(b,2)$in( + )$paren(5$in( )sin$paren(5$in( )%tau)$in( - )132$in( )sin$paren(3$in( )%tau)$in( + )599$in( )sin$paren(%tau),$(,$))$in( )$sup(a,4)$in( )b$in( + )$paren(cos$paren(5$in( )%tau)$in( - )24$in( )cos$paren(3$in( )%tau)$in( + )23$in( )cos$paren(%tau),$(,$))$in( )$sup(a,5),$(,$))$in( )$sup(e,2),1024)$hinge()$in( - )$q($paren($paren(sin$paren(3$in( )%tau)$in( + )9$in( )sin$paren(%tau),$(,$))$in( )$sup(b,3)$in( + )$paren(3$in( )cos$paren(3$in( )%tau)$in( - )3$in( )cos$paren(%tau),$(,$))$in( )a$in( )$sup(b,2)$in( + )$paren(21$in( )sin$paren(%tau)$in( - )3$in( )sin$paren(3$in( )%tau),$(,$))$in( )$sup(a,2)$in( )b$in( + )$paren(cos$paren(%tau)$in( - )cos$paren(3$in( )%tau),$(,$))$in( )$sup(a,3),$(,$))$in( )e,32)$hinge()$in( + )sin$paren(%tau)$in( )b$hinge()$in( + )cos$paren(%tau)$in( )a$close(])$ina($, )$hinge()%tau$hinge()$in( = )$open($()$in( - )$q($paren(69$in( )$sup(b,4)$in( + )138$in( )$sup(a,2)$in( )$sup(b,2)$in( + )21$in( )$sup(a,4),$(,$))$in( )$sup(e,2),256)$hinge()$in( + )$q($paren(3$in( )$sup(b,2)$in( + )3$in( )$sup(a,2),$(,$))$in( )e,8)$hinge()$in( + )1$close($))$hinge()$in( )t$close(])$close(])
  110.  
  111. (c10) /* 
  112.    Finally, consider the Duffing oscillator with weak sinusoidal forcing. */
  113.  
  114. forced_duffing: 'diff(x,t,2)+x+e*x^3=e*gamma*cos(t);
  115. |$label(0,15,Times New Roman,$(d10$))$q($sup(d,2)x,d$sup(t,2))$hinge()$in( + )e$in( )$sup(x,3)$hinge()$in( + )x$hinge()$in( = )e$hinge()$in( )cos$paren(t)$hinge()$in( )$greektext(g)
  116.  
  117. (c11) lindstedt(forced_duffing,e,1);
  118. |$label(0,15,Times New Roman,$(d11$))$open([)$open([)$open([)$q($paren(3$in( )cos$paren(3$in( )%tau)$in( - )cos$paren(%tau),$(,$))$in( )e$in( )$greektext(g),72)$hinge()$in( + )$q($sup(3,2$in(/)3)$in( )$sup(4,1$in(/)3)$in( )cos$paren(%tau)$in( )$sup($greektext(g),1$in(/)3),3)$close(])$ina($, )$hinge()%tau$hinge()$in( = )t$close(])$close(])
  119.